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Abstract 

In synthetic aperture radar (SAR), images are formed by focusing the response of stationary objects to a single 
spatial location. On the other hand, moving targets cause phase errors in the standard formation of SAR images 
that cause displacement and defocusing effects. SAR imagery also contains significant sources of non- stationary 
spatially-varying noises, including antenna gain discrepancies, angular scintillation (glints) and complex speckle. In 
order to account for this intricate phenomenology, this work combines the knowledge of the physical, kinematic, and 
statistical properties of SAR imaging into a single unified Bayesian structure that simultaneously (a) estimates the 
nuisance parameters such as clutter distributions and antenna miscalibrations and (b) estimates the target signature 
required for detection/inference of the target state. Moreover, we provide a Monte Carlo estimate of the posterior 
distribution for the target state and nuisance parameters that infers the parameters of the model directly from the 
data, largely eliminating tuning of algorithm parameters. We demonstrate that our algorithm competes at least as 
well on a synthetic dataset as state-of-the-art algorithms for estimating sparse signals. Finally, performance analysis 
on a measured dataset demonstrates that the proposed algorithm is robust at detecting/estimating targets over a 
wide area and performs at least as well as popular algorithms for SAR moving target detection. 

I. Introduction 

This work proposes algorithms for detecting and estimating targets in synthetic aperture radar (SAR) 
images. The image formation process for SAR images is more complicated than that of standard electro- 
optical images. Examples of these complexities include: 

• SAR images have complex- valued rather than real-valued intensities, and the SAR phase information 
is of great importance for detection and estimation of target states. [l]-[3]. 

• SAR images are corrupted by spatiotemporally-varying antenna gain/phase patterns that often need 
to be estimated from homogeneous target-free data [4], [5]. 

• SAR images have spatially-varying clutter that can mask the target signature unless known a priori 
or properly estimated [6]. 

• SAR images contain motion-induced displacement and diffusion of the target response [1], [7]. 

• SAR images include multiple error sources due to radar collection and physical properties of the 
reflectors, such as angular scintillation (a.k.a. glints) [8] and speckle [9], [10]. 

Despite these complications, a great deal of structure exists in SAR images that can be leveraged to 
provide stronger SAR detection and tracking performance. This includes (a) using the coherence between 
multiple channels of an along-track radar in order to remove the stationary background (a.k.a, 'clutter'), 
(b) assuming that pixels within the image can be described by one (or a mixture) of a small number of 
object classes (e.g., buildings, vegetation, etc.), and (c) considering kinematic models for the target motion 
such as Markov smoothness priors. From this structure in SAR imagery, one might consider models that 
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assume that the clutter Ues in a low-dimensional subspace that can be estimated directly from the data. 

Indeed, recent work Borcea et al. [11] has shown that SAR signals can be represented as a composition 
of a low-rank component containing the clutter, a sparse component containing the target signatures, and 
additive noise. 

In general, SAR images are formed by focusing the response of stationary objects to a single spatial 
location. Moving targets, however, will cause phase errors in the standard formation of SAR images that 
cause displacement and defocusing effects. Most methods designed to detect the target depend on either 

(a) exploiting the phase errors induced by the SAR image formation process for a single phase center 
system or (b) canceling the clutter background using a multiple phase center system. In this chapter, we 
provide a rich model that can combine (and exploit) both sources of information in order to improve on 
both methodologies. 

Fienup [1] provides an analysis of SAR phase errors induced by translational motions for single-look 
SAR imagery. He shows that the major concerns are (a) azimuth translation errors from range-velocities, 

(b) azimuth smearing errors due to accelerations in range, and (c) azimuth smearing due to velocities 
in azimuth. Fienup also provides an algorithm for detecting targets by their induced phase errors. The 
algorithm is based on estimating the moving target's phase error, applying a focusing filter, and evaluating 
the sharpness ratio as a detection statistic. Jao [7] shows that given both the radar trajectory and the target 
trajectory, it is possible to geometrically determine the location of the target signature in a reconstructed 
SAR image. Although the radar trajectory is usually known with some accuracy, the target trajectory is 
unknown. On the other hand, if the target is assumed to have no accelerations, Jao provides an efficient 
FFT-based method for refocusing a SAR image over a selection of range velocities. Khwaja and Ma 
[12] provide a algorithm to exploit the sparsity of moving targets within SAR imagery; they propose a 
basis that is constructed from trajectories formed from all possible combinations of a set of velocities and 
positions. To combat the computational complexity of searching through this dictionary, the authors use 
compressed sensing techniques. Instead of searching over a dictionary of velocities, our work proposes 
to use a prior distribution on the target trajectory that can be provided a priori through road and traffic 
models or adaptively through observations of the scene over time. 

The process of removing the stationary background in order to detect moving targets is also known in 
the literature as 'change detection' or 'clutter suppresion'. GieruU [13] provides a statistical analysis of 
the phase and magnitude of complex SAR images for two channels. He shows that SAR images cannot be 
modeled as spatially-invariant Gaussian in many cases of interest, such as in urban environments, where 
the statistics vary spatially and may be modulated by random variations. In our work, we model the 
distributions of the clutter as spatially varying and model the random modulations directly. 

Ender [6] applies space-time adaptive processing (STAR) to multiple-channel SAR imagery. Similar 
to standard change detection algorithms such as displaced phase center array (DPCA) and along-track 
interferometry (ATI), STAR models the clutter as being embedded in a one-dimensional subspace. However, 
STAP extends those algorithms to using N > 2 channels, where a single channel is used to estimate the 
stationary background and the remaining (N — 1) channels are used to estimate the moving component. 
However, STAP relies on estimating the complex-valued covariance matrix of the A^^-channel system, 
which in turn depends on the availability of homogeneous target-free secondary data. 

There are a multitude of algorithms for change detection that are based on multi-temporal SAR images 
rather than multi-channel data. Bazi and Bruzzone [14] develop methods for multi-temporal change 
detection that use adaptive thresholds for declaring changes based on a theoretical analysis of a generalized 
Gaussian model. Bovolo and Bruzonne [15] provide another algorithm for change detection that employs 
a wavelet-based multiple scale decomposition of multitemporal SAR images, with an adaptive scale driven 
fusion algorithm. 

Ranney and Soumekh [4], [5] develop methods for change detection from SAR images collected at 
two distinct times that are robust to errors in the SAR imaging process. They address error sources 
including inaccurate position information, varying antenna gains, and autofocus errors. They propose 
that the stationary components of multi-temporal SAR images can be related by a spatially-varying 2- 
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dimensional filter. To make the change detection algorithm numerically practical, the authors propose that 

this filter can be well- approximated by a spatially invariant response within small subregions about any 
pixel in the image. This thesis adopts this model for the case where there are no registration errors. Under 
a Gaussian assumption for the measurement errors, it can be shown that the maximum likelihood estimate 
for the filter coefficients can be computed easily through simple least squares. 

Ground Moving Target Indication (GMTI) methods involve the processing of SAR imagery to detect 
and estimate moving targets. Often clutter cancellation and change detection play a preprocessing role in 
these algorithms [16]-[19]. This chapter aims to combine properties of many of these algorithms into a 
unifying framework that simultaneously estimates the target signature and the nuisance parameters, such 
as clutter distributions and antenna calibrations. 

It should be noted that many of the previously discussed algorithms work well in certain situations, 
but do not provide estimates of their uncertainty that may be necessary for adaptive sensing, sensor 
management, or sensor fusion. This chapter aims to bridge this gap by providing a Bayesian formulation 
that provides uncertainty distributions for the presence of the moving targets and their positions. Under this 
Bayesian formulation, we can generate the posterior distribution of the target state(s) given the observations 
(i.e., the SAR images). 

Recently, there has been great interest by Wright et al. [20], Lin et al. [21], Candes et al. [22] and 
Ding et al. [23] in the so-called robust principal component analysis (RPCA) problem that decomposes 
high-dimensional signals as 

I^L + S + E, (1) 

where / G M^^^ is an observed high dimensional signal, L G ]R^^^ is a low-rank matrix with rank 
r < NM, S G M^x^ is a sparse component, and E G M^^^ is dense low-amplitude noise. In [20]-[22], 
inference in this model is done by optimizing a cost function of the form 

argmin + 7 ||5||, + ||/ - - S\\^ (2) 

where the last term is sometimes replaced by the constraint I — L + S. One major drawback of these 
methods involves finding the algorithm parameters (e.g., tolerance levels or choices of 7, /i), which may 
depend on the given signal. Moreover, it has been demonstrated that the performance of these algorithms 
can depend strongly on these parameters. 

Bayesian methods by Ding et al. [23] have been proposed that simultaneously learn the noise statistics 
and infer the low-rank and sparse components. Moreover, they show that their method can be generalized 
to richer models, e.g. Markov dependencies on the target locations. Additionally, these Bayesian inferences 
provide a characterization of the uncertainty of the outputs through a Markov Chain Monte Carlo (MCMC) 
estimate of the posterior distribution. 

The work by Ding et al. [23] is based on a general Bayesian framework [24] by Tipping for obtaining 
sparse solutions to regression and classification problems. Tipping's framework uses simple distributions 
(e.g., those belonging to the exponential class) that can be described by few parameters, known as 
hyperparameters. Moreover, Tipping considers a hierarchy where the hyperparameters themselves are 
assumed to have a known 'hyperprior' distribution. Often the prior and hyperprior distributions are chosen 
to be conjugate, so that inference is simple. Tipping provides insight into choosing the hyperparameter 
distributions so as to be non-informative with respect to the prior. This latter property is important in 
making it possible to implement inference algorithms with few tuning parameters. Finally, Tipping provides 
a specialization to the 'relevance vector machine' (RVM), which can be thought of as a Bayesian version 
of the support vector machine. Wipf et al. [25] provides an interpretation of the RVM as the application of 
a variational approximation to estimating the true posterior distribution. Wipf et al. explains the sparsity 
properties of the sparse Bayesian learning algorithms in a rigorous manner. Additionally, it also provides 
connections with other popular work in sparse problems, such as the FOCUSS and basis pursuit algorithms. 

We adopt this hierarchical Bayesian model to SAR images. This requires the following non-trivial 
extensions: (a) we consider complex-valued data rather than real-valued intensity images; (b) we model 
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correlated noise sources based on physical knowledge of SAR phase history collection and image forma- 
tion; (c) we relax the assumption of a low-rank background component by assuming that the background 
component lies in a low-dimensional subspace; and (d), we directly model SAR phenomenology by 
including terms for glints, speckle contributions, antenna gain patterns, and target kinematics. Moreover, 
we demonstrate the performance of the proposed algorithm on both simulated and measured datasets, 
showing competitive or better performance in a variety of situations. 

The rest of the paper is organized as follows: Notation is given in Section II and the image model 
is provided in III. Markov, spatial, and/or target kinematic extensions are discussed in Section IV. The 
inference algorithm is given in Section V. Performance is analyzed over both simulated and measured 
datasets in Section VI. We conclude and point to future work in Section VII. 

II. Notation 

Available is a set of SAR images of a region formed from multiple passes of an along-track radar 
platform with multiple antennas (i.e., phase centers.) Moreover, images are formed over distinct azimuth 
angle ranges that can be indexed by the frame number, /. Table I provides the indexing scheme used 
throughout this chapter in order to distinguish between images from various antennas, frames, and/or 
passes. Table 11 provides a list of indexing conventions used to denote collections of variables. 

TABLE I 

Index variable names used in paper 



Index Description 


Index Variable 


Range 


Antenna (channel) 


k 


1,2,. 


.,K 


Frame (azimuth range) 


f 


1,2,. 


■,F 


Pass 


i 


1,2,. 


.,N 


Pixel 


P 


1,2,. 


■,P 



TABLE n 
Our data indexing conventions 



Variable 


Convention 


Description 


Ap) 


Standard 


Value at pixel p, antenna k, 

and frame /, pass i 


Ap) 


Underline 


Values at pixel p, frame /, 
and pass i over all antennas 


Ap) 


Lower-case, 
Boldface 


Values at pixel p and frame / 
over all antennas and passes 




Upper-case 
Boldface 


Values over all pixels and 
antennas at frame / and pass i 


I 


Upper-case, 
Boldface, No Indices 


Values over all pixels, antennas, 
frames, and passes 



We model the quadrature components of the SAR images with the complex-normal distribution, where 
we use the notation 

w^cu (0, r) (3) 

where CN'{fJ,, T) represents the complex-Normal distribution with mean n and complex covariance matrix 
r, and w is random vector of K complex-values (from each of K antennas.) 
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Fig. 1 . This figure provides a graphical representation of the proposed SAR image model. The dark circle represents the observed random 
variable. The unshaded circles represent the basic parameters of the model, while the dashed circles represent hyperparameters that are also 
modeled as random variables. 

III. SAR IMAGE MODEL 

We propose a decomposition of SAR images at each frame / and pass i as follows 

If,, = Hf,, o {Lf, + Sf, + Vf,,) , (4) 

where Hf,i is a spatiotemporally-varying filter that accounts for antenna calibration errors, L f,i is a low- 
dimensional representation of the background clutter, j is a sparse component that contains the targets 
of interest, Vf,i is zero-mean additive noise, and o denotes the Hadamard (element-wise) product. Each 
of these components belongs to the space C^^^. The remainder of this section discusses the model in 
detail. Figure HI shows a graphical representation of the model. 

A. Low-dimensional component, Lf,, 

We propose a decomposition of the low-rank component as 

where Bf is the inherent background that is identical over all passes, X/ j is the speckle noise component 
that arises from coherent imaging in SAR. Posner [9] and Raney [10] describe speckle noise, which tends 
to be spatially correlated depending on the texture of the surrounding pixels. 

The quadrature components of radar channels are often modeled as zero-mean Gaussian processes, 
though GieruU [13] demonstrates that for heterogeneous clutter (such as in urban scenes), one must 
consider spatially varying models. To account for this spatial variation, this model assumes that each 
background pixel can be defined by one of J classes that may be representative of roads, vegetation, or 
buildings within the scene. Our model is low-dimensional since J <^ P, where P is the number of pixels 
in the measured images. We put a multinomial model on each object class 

cip) ^ [cf]'^^ ~ Multinomial(l;?i,?2,...,gj) (6) 

where qj is the prior probability of the j-th object class. Then the class assignment C^^^ is the single 
location in c with value equal to one. We use a hidden Markov model dependency that reflects that 
neighboring pixels are likely to have the same class. The class C^^'^ defines the distribution of the pixel 
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p, where we specifically model the background and speckle components respectively as complex-normal 
distributed: 

bf - CAA (o, rg'"^ ) , xfj^ CN (o, r?"^ ) (7) 

Note that the class type specifies the distribution of the pixels and each vector of K values (e.g. background 
0'' or speckle x^']) is drawn independently from that distribution. 

B. Sparse component, Sf^i 

The sparse component is modeled as 

Sf, = (A^ ® ll) o Gj, + (Ajf, ® ll) o Mf,, (8) 

where Gf^i G C^^^ is the specular noise (glints) component with associated indicator variables e 
{0, 1}^, Mf.i G C^^^ is the (moving) target component with associated indicator variables A^- G 

{0, 1}^, Ik is the all ones vector of size K x 1, and ® is the Kronecker product. Note that this shared 
sparsity model assumes that the glint/target components are present in one antenna if and only if they 
are present in the other antennas. Moreover, glints are known to have a large angular dependence, in the 
sense that the intensity of the glint dominates in only a few azimuth anglesbut is present from pass to 
pass as described by Borden [8]. Thus, the indicators for glints do not depend on the pass index i. Once 
again, we assume that the glints and target components are zero-mean complex-normal distributed with 
covariances Tq and Tm, respectively. 

The indicator variable 5^'^'^^ at pixel p where z is representative of either or m is modeled as 

S''^^ ~ Bernoulli(7r^'(P)), (9) 
7r^'(f) ~ Beta(a^, b^) (10) 

A sparseness prior is obtained by setting a-^/la-j^ + feTr] <^ 1- Alternatively, we can introduce additional 
structure in our model by letting and bj, depend on previous frames (temporally) and/or neighboring 
pixels (spatially). This is particularly useful for detecting multi-pixel targets that move smoothly through 
a scene. Section IV discusses this modification in greater detail. 



C. Distribution of quadrature components 

Many SAR detection algorithms rely on the ability to separate the target from the background clutter 
by assuming that the clutter lies in a low-dimensional subspace of the data. Consider a random vector 
of complex variables w ~ CJ\f (0, F) where w is representative of b, x, g or m. Under the assumptions 
that (a) the quadrature components of each antenna are zero-mean normal with variance and (b) the 

then r can be shown to have the form 

-j<t>12 . . . /-)f>~^<t>lK 



correlation among components Wm and Wn is given by pe 



pe 



1 

J</>12 



pe 



pe 
pe 



-j<t>2K 



.j<t>lK 



pe- 



.j<t>2K 



(11) 



where cr^ is the channel variance, p is the coherence between antennas, and {0nm}„ ^ are the inter- 

ferometric phase differences between the antennas ^ In an idealized model with a single point target, 
the interferometric phases can be shown to be proportional to the target radial velocity. In images 
containing only stationary targets (i.e., the background components), the covariance matrix has a simpler 
form: 

'^backgrmmd — ((1 " P)IkxK + P'^K^Jc) (12) 



'a more general model could account for different channel variance and coherence values, but since we use the calibration constants Jf/,, 
to equalize the channels, the effect was seen to be relatively insignificant. 
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where Irxk is the X x X identity matrix and Ik is the all-ones vector of length K. 

It should be noted that the covariance matrix in equation (11) is related directly to some common 
methods for change detection in SAR imagery. In particular, consider the two antenna case (K = 2). 
Along-track interferometry (ATI) thresholds the phase (pu in order to detect moving targets which have 
non-zero phases. Moreover, one can easily show that the eigendecomposition of T leads to eigenvalues A 
and eigenvectors u: 

A(r) = {2(7^(1 + 2(7^(1 -p)} (13) 
^ ^ ■ (14) 



r 



1 



Displaced phase center array (DPCA) processing thresholds the difference between the two channels. 
Indeed, for small phases, the second eigenvector of T reduces to [1; —1]^. Thus DPCA can be interpreted 
as a projection onto the eigenvector of T. Deming [2] shows that ATI performs well when canceling bright 
clutter (i.e., high cr^ and p fti 1), while DPCA performs well for canceling dim clutter (i.e., small cr^ and 
p ^ 0.) In our work, we combine the discriminating power of both DPCA and ATI by modeling the 
covariance matrices directly. Ender [6] provides space-time adaptive processing (STAP), where optimal 
detection schemes for moving targets are based on the estimation of T. However, the performance of STAP 
depends on the availability of target-free homogeneously distributed measurements in order to estimate 
r effectively. In this chapter, we simultaneously estimate the covariance matrices as well as the target 
contributions. Thus, we demonstrate the capability to detect targets even in the presence of heterogeneous 
measurements. 

In this thesis, the covariance matrix T is modeled as a random variable using a modified version of the 
Multivariate-Normal-Inverse-Wishart conjugate distributions. In particular, we let 



w 



CAr{o,(T^rp) (15) 

Tp ~ InvWishart (ar((l - p)Ikxk + pIrIr), i^r) (16) 
(j^ ~ InvGamma(ao-, 6(t) (17) 
p ~ Bcta(ap, bp) (18) 

where a„ = h„ = 10^^ as suggested by Tipping [24] to promote non-informative priors, (a^, hp) are chosen 
so that p ~ 1 to ensure a high coherence among the background components, z/p is a parameter that controls 
how strongly to weight the prior covariance matrix, and or is chosen so that E\r = (1— p)/xx x+plft:!]^- 
In this thesis, we choose vy to be large to reflect our belief that a'^Tp should be close to equation (11). Note 
that this model separates the learning of the channel variance cr^, which we have no a priori knowledge 
about, from the learning of the correlation structure Tp. 

D. Calibration filter, Hf^i 

The caUbration constants are assumed to be constant within small spatial regions p G Zg, though they 
may vary as a function of antenna, frame, or pass. In particular, we let 

^S,. = ^fcj.(^?),VpeZ„ (19) 

ZkjM-CM{l,{a")^) (20) 

where we note that if (cr^)^ is large, then maximum likelihood inference in this case yields the least- 
squares solution. 
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E. Summary of SAR Image Model 

Tables HI provides a summary of the distributions for the proposed decomposition of SAR images. The 

table also provides a characterization of spatial (across pixels) and temporal (across frames and passes) 
dependencies. For example, background and speckle components have distributions characterized by their 
class j. Thus, all pixels with class j belong to a subset Qj C {1, 2, ... , P}. In contrast, the distribution of 
moving targets is assumed to be identical across all pixels, yet the distribution of their indicators varies 
for each pixel, frame, and pass. 

Tables IV and V provide a summary of the parameters of the distributions in Table III. We provide 
the simple model for target and glint indicator probabilities that just assumes that they are sparse in the 
image. We can introduce additional richness in the model by allowing the parameters a^^ and 6^ to vary 
over pixels, frames, and passes as described in Section IV-A. 



IV. Markov/spatial/kinematic models for the sparse component 
A. Indicator probability models 

This model contains multiple indicator variables with prior probabilities distributed as Bcta(a7r, ^vr)- 
Moreover, sparsity is obtained when a^/[a7r + 67r] <^ 1. Alternatively, we can introduce additional structure 
in our model by letting a-j^ and bj^ depend on previous frames (temporally) and/or neighboring pixels 
(spatially). This is especially useful for detecting multi-pixel targets that move smoothly through a scene. 

Define {p, A^-) to be a function that maps the indicator variables A^,- to a real number. For 
example, this may be the average number of non-zero indicators in the neighborhood of pixel p, or a 
weighted version that puts higher value on neighboring pixels. For / = 1, we let 



and for / > 1 



KAp) 




W^{p,Af;^)>e^^^,,^„ 
else, 



(21) 



[aH hnf, W^ip 



[ttL hi 



W^{p, Af 
else. 



> e 



M 

spatial 



and 



) > 

1,1/ temporal 1 



(22) 



In this chapter, we choose {aL,bL, aH,bH) so that ai/iai ^-b + L) <C 1 and an/icLH + h + H) >> 0. A 
similar model can be introduced for the probabilities of the glints. 



B. Target kinematic model 

In some apphcations, such as target tracking or sequential detection, we may have access to an estimate 
of the kinematic state of the target(s) of interest, such as position, velocity and acceleration. This may be 
useful for predicting the location of the target at sequential frames. For simplicity, consider a single target 
at time r whose state ^(r) = {r{T),r{T)) is known with standard errors S^(r). Note that the uncertainty 
model for (r, r) may be (a) known a prior from road maps or traffic behavior patterns, or (b) learned 
adaptively using some signal processing algorithm such as the Kahnan or particle filters. 

In standard SAR image formation, moving targets tend to appear displaced and defocused as described 
by Fienup [1] and Jao [7]. Moreover, Jao showed that given the radar trajectory {q.q) and the target 
trajectory {r,r), one can predict the location of the target signature within the image p by solving a 
system of equations that equate Doppler shifts and ranges, respectively, at each pulse: 

^[\\p-q{r)h-\\r{r)-q{r)U^^^.=Q (23) 
\\p*-q{T)\\,^\\r{T)-q{T)\\,, (24) 



Temporal 
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All /, All i 


All /, All i 


All /, All i 
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procedure {0} 


= SARGibbs(0o, J) 

A'-i * samples 










for iteration 


1 to N ^urnin ~l~ samples dO 




Sample ~ 


/(S,X,G,M,A«,A^|J,-) 


//Base 


Sample ~ 


f{H\I,-) 


//Calibration filter 


Sample ~ 


/(C|/,-) 


//Class assignment 


Sample ~ 


J yn\^^ -) 


II Hyper-parameters 


® iteration- 


Nburnin ^ ^ iteTatiOfl > Nburnin 




end for 






end procedure 







Fig. 2. Gibbs Sampling Pseudocode 



which can be reduced to the simpler system of equations: 

q{T) ■ [p* - q{T)] = [r(r) - g(r)] ■ [r(r) - q(r)] (25) 
\\p*-q(r)\\, = \\r{r)-q{r)\\, (26) 

The probable locations of the target can be predicted by one of several methods, including: 

• Monte Carlo estimation of the target posterior density. 

• Gaussian approximation using linearization or the unscented transformation to approximate the pos- 
terior density 

• Analytical approximation. 

Given an estimate of the posterior density, we can modify the function described in the previous 
section to include dependence on this kinematic information. Details of the posterior density estimation 
are provided in Appendix A. 

V. Inference 

In the proposed hierarchical model, the distribution of hyper-parameters at the base layer are generally 
chosen to be conjugate to the distributions at the next layer. This allows for efficient approximation methods 
for the posterior distribution in the sense that we can sample exactly from these distributions. In particular, 
we use a Markov Chain Monte Carlo (MCMC) algorithm in the form of a Gibbs sampler to iteratively 
estimate the full joint posterior. In MCMC, this distribution is approximated by drawing samples iteratively 
from the conditional distribution of each (random) model variable given the most recent estimate of the 
rest of the variables (which we denote by — ). Let © = {B,X, G, M, A*^, A^^, H, C, ry} represent a 
current estimate of all of the model variables where r] represents the set of all hyper-parameters. Given 
measurements /, the inference algorithm is given in Figure 2. Note that MCMC algorithms require a 
burn-in period after the Markov chain has become stable, where the duration of burn-in period depends 
on the problem. After this point, we collect Ngampies samples that represent the full joint distribution. 
However, we point out a couple of important features here. First, the sampling of the base model can be 
rewritten as 

/(S,X,G,M,A^,A^^|/,-) (27) 

J \y.f )"^/,l:iV'y/,l:Ar)''*'/,l;-'V' / ' I ' ) 

P,f 

The conditional independence among pixels and frames given the nuisance parameters allows us to easily 
parallelize the sampling procedure over the largest dimensions of the state. Moreover, we can extend the 
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parallelization to sampling independently over passes by separating the sampling of equation (27) into 
two Gibbs steps from the densities: 

/(6,^\x,^L^,m,^)^,<5,^;|^)|J,-) (28) 

i 

f{9%.N,5f^'\l,-) (29) 

i 

In both of the sampling steps in equations (28) and (29), we have an exact inference algorithm over 
multivariate-Gaussian distributed variables and Bernoulli distributed variables. This leads to faster con- 
vergence of the Markov chain and subsequently fewer bum-in samples. The conditional density for the 
nuisance parameters rf given the remainder variables can also be re-written to allow for efficient sampling. 
In particular, due to conditional independence we have: 

/(r/l J, -) =/(r^|M, A^)/(7r^|M, A^, V^) 

•/(r«|G,A«)/(7r«|G,A«,r«) ^3^^ 

■\[f{T^mB.C)f{T^{j)\X,C) 

3 

where T represents the parameters related to the covariance matrices (i.e., the variance u^, correlation 
structure T p, and the coherence p). Once again, this decomposition allows for a sampling procedure that 
leads to faster convergence of the Gibbs sampler. Moreover, the sampling procedures for the individual 
densities in equation (30) tend to require sufficient statistics that are of significantly smaller dimension 
and thus more desirable from a computational viewpoint. For example, sampling of the covariance matrix 
depends only on a X x sample covariance matrix. It should be noted that sampling of the 
covariance matrices requires additional effort in order to constrain its shape to that of equation (11). 
In particular, we use a Metropolis-Hastings step, which can be easily done by noting that the posterior 
density /(r^,p^, ((T^)^|Ty) is proportional to an Inverse-Wishart distribution. Details are provided in 
Appendix B. 



VI. Performance analysis 

A. Simulation 

We first demonstrate the performance of the proposed algorithm, which we refer to as the Bayes SAR 
algorithm, on a simulated dataset. Images were created according to the model given in Section HI with 
parameters given in Table VI. The low-dimensional component was divided into one of two classes ('dim' 
or 'bright'). Pixels were deterministically assigned to one of these classes to resemble a natural SAR image 
(see Figure 4). The sparse component included a randomly placed target with multiple-pixel extent. A 
spatiotemporally varying antenna gain filter was uniformly drawn at random on the range [0, 27r) for 
groups of pixels of size 25 x 25. Lastly, zero-mean IID noise was added with variance Onoise- 

The Bayes SAR model is applied to infer the low-dimensional component Lf^i and sparse target 
component S/^ with estimates denoted Lf^ and Sj^i, respectively. Hyperparameters of the model are 
chosen according to the Section V. Results are given by the mean of MCMC inference with 500 bum-in 

iterations followed by 100 collection samples. We consider three metrics to evaluate the reconstruction 

\\s-s\\ \\s-s\\ , , . , , . , 

errors: ^ , u„u ^ , u„u ^ , where the norm is taken over the vectorized quantities. 

Il-'^lb ll'-'lh IPIIo 

In comparison to the Bayes SAR model, results are given for state-of-the-art algorithms for Robust 
Principal Component Analysis (RCPA): an optimization-based approach proposed by Wright et al. [20] 
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TABLE VI 
Parameters of simulated dataset 



Parameter 


Value 


Pixels in image, P 


P = 100 X 100 


Number of frames per pass, F 


F = l 


# of antennas, K 


K = 3 


# of passes, A'" 


N e {5,10,20} 


# of target pixels/image, Ntargets 


Ntargets = 20 


Clutter of background, p 


p€ {0.9,0.99,0.999,0.9999} 


Variance of targets, a^arget 


2 _ 1 

^target 


Variance of background 


Either aji^ = a^,„„<;^/100 

2 2 
or CTbright — '^clutter 


Signal-to-noise-plus clutter (SCNR) 


SCNR = , 

clutter ' noise 

e {0.1,0.5,1,2} 



TABLE Vll 

Comparison of proposed method (Bayes SAR) to RPCA Methods with = 20, F = 1, 7^ = 3. Note that the Bayes SAR 
method performs about twice as well as either of the rpca methods for all criteria. the bayes sar method also 

produces a sparse result. 



(a) Bayes SAR 



(b) Opt. RPCA 



SCNR 


Coherence 








I|i|l2 


l|S|l2 


IISIIo 


10% 


0.900 


0.057 


0.578 


0.550 


10% 


0.9999 


0.045 


0.419 


0.367 


100% 


0.900 


0.057 


0.155 


0.150 


100% 


0.9999 


0.052 


0.122 


0.096 


200% 


0.900 


0.057 


0.123 


0.137 


200% 


0.9999 


0.056 


0.114 


0.092 



(c) Bayes RPCA 



SCNR 


Coherence 


ll^-^ll. 






\\L\\2 


11S112 


IISIIo 


10% 


0.900 


0.117 


0.998 


3.761 


10% 


0.9999 


0.108 


0.990 


3.799 


100% 


0.900 


0.116 


0.764 


3.451 


100% 


0.9999 


0.117 


0.741 


3.494 


200% 


0.900 


0.125 


0.706 


3.665 


200% 


0.9999 


0.129 


0.692 


3.720 



SCNR 


Coherence 








\\L\\2 


11SII2 


IISIIo 


10% 


0.900 


0.111 


3.175 


111.026 


10% 


0.9999 


0.113 


3.237 


109.716 


100% 


0.900 


0.111 


1.189 


109.520 


100% 


0.9999 


0.110 


1.173 


108.203 


200% 


0.900 


0.112 


1.058 


111.120 


200% 


0.9999 


0.110 


1.035 


109.583 



and Candes et al. [22] and a Bayesian-based approach proposed by Ding et al. [23]^. The optimization- 
based approach requires a tolerance parameter which is related to the noise level, as suggested by Ding 
et al. [23]. We chose this parameter in order to have the smallest reconstruction errors. The Bayesian 
method did not require tuning parameters, except for choosing the maximum rank of L f^i which was set 
to 20. 

||5_5|| 

Figure 3 compares the relative reconstruction error of the sparse (target) component, u^.. \ across 
all algorithms, number of passes N, coherence of antennas p, and SCNR. In all cases, the Bayes SAR 
method outperforms the RPCA algorithms with improving performance if either coherence or SCNR 

^For the optimization-based approach, we used the exact_alm_rpca package (MATLAB) by Lin et al. [21], downloaded from http://watt. 
csl.illinois.edu/perceive/matrix-rank/home.html. For the Bayesian-based approach, we used the Bayesian robust PCA package, downloaded 
from http://www.ece.duke.edu/~lihan/brpca_code/BRPCA.zip. 
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Fig. 3. This figure compares thie relative reconstruction error of the target component, ^ , as a function of algorithm, number of passes 

N, coherence of antennas p, and signal-to-clutter-plus-noise ratio (SCNR). From top-to-bottom, the rows contains the output of the Bayes 
SAR algorithm, the optimization-based RPCA algorithm, and the Bayes RPCA algorithm. From left-to-right, the columns show the output 
for N = 5, N = 10, and A'^ = 20 passes (with F = 1 frames per pass). The output is given by the median error over 20 trials on a 
simulated dataset. It is seen that in all cases, the Bayes SAR method outperforms the RPCA algorithms. Moreover, the Bayes SAR algorithm 
performs better if either coherence increases (i.e., better clutter cancellation) or the SCNR increases. On the other hand, the performance of 
the RPCA algorithms does not improve with increased coherence, since these algorithms do not directly model this relationship. 



increases. Table VII provides additional numerical results for the case = 20. The RCPA algorithms 
perform poorly in reconstructing the sparse component with relative errors near or greater than 1. This 
reflects the fact that (a) these algorithms miss significant sources of information, such as the correlations 
among antennas and among quadrature components, and (b) = 20 may be too few samples to reliably 
estimate the principal components in these non-parametric models. In measured SAR imagery, it might 
be unreasonable to expect ^ 20 passes of the radar, which suggests that these RPCA algorithms 
will likely perform poorly on such signals. In contrast, it is seen that the Bayes SAR method obtains 
low reconstruction errors for both low-dimensional and sparse components as either coherence or SCNR 
increase. 
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(a) L + S (b) L (c) S 



Fig. 4. This figure provides a sample image used in the simulated dataset for comparisons to RPCA methods, as well as its decomposition 
into low-dimensional background and sparse target components. This low SCNR image is typical of measured SAR images. Note that the 
target is randomly placed within the image for each of A'^ passes. In some of these passes, the target is placed over low-amplitude clutter 
and can be easily detected. In other passes, the target is placed over high-amplitude clutter, which reduces the capability to detect the target. 



B. Measured data 

In this section, we compare performance of the Bayes SAR approach using a set of measured data from 
the 2006 Gotcha SAR sensor collection. In particular, images were formed from phase histories collected 
over a scene of size 375m by 1200m for = 3 passes and K = 3 antennas. Each image was created 
with a coherent processing time of 0.5 seconds with the addition of a Blackman-Harris window in the 
azimuth direction to reduce sidelobes. Images were created at overlapping intervals spaced 0.25 seconds 
apart for a total of 18 seconds. Note that the ability to take advantage of correlated images (as in this 
case) is one of the benefits of using the proposed model/inference algorithm. 

We consider three alternative approaches in comparison to the Bayes SAR approach: (1) displaced-phase 
center array (DPCA) processing, (2) along-track interferometry (ATI), and (3) a mixture of DPCA/ATI. 
Note that all variants of ATI/DPCA depend on the chosen thresholds for phase/magnitude, respectively. 

1) Comparisons to DPCA/ATI: We begin by comparing the output of the proposed algorithm across the 
entire 375m by 1200m scene. Figure 5 shows the output of the Bayes SAR algorithm and the DPCA/ATI 
comparisons. It is seen that there are significant performance gains by using calibrated images as shown 
in (c) and (f) as compared to their original versions, (b) and (e), respectively. Furthermore, the proposed 
approach also provides a sparse output without choosing thresholds as required by DPCA/ATI. Note that 
in this figure, calibration is accomplished by using the outputs Hf^i from the Bayes SAR approach. 

Figures 6 and 7 display the detection performance over two smaller scenes of size 125m by 125m 
as a function of magnitude and phase, respectively. For each scene, images are provided for sequential 
scenes separated by 0.5 seconds. Scene 1 contains strong clutter in the upper left region, while Scene 2 
has relatively little clutter. It is seen that the proposed approach (2nd column) provides a sparse solution 
containing the targets of interest in each of the 4 images. Moreover, the 3rd column provides the estimated 
probability that a target occupies a given pixel, in comparison to the (0,1) output of DPCA. Although 
most estimated probabilities are near 1, there are a few cases where this is not the situation: in scene 
2(d), a low-magnitude target is detected with low probability in the lower-right; in scene 1(b) a few target 
pixels from the clutter region are detected with low probability. In contrast, the performance of DPCA 
depends strongly on the threshold. In Scene 1, a 30 dB threshold provides a large number of false alarms. 
However, in Scene 2, the low-magnitude targets are missed for the 15 dB threshold, but detected at the 
30 dB threshold. 

Figure 7 shows the detection performance based on phase over the same 4 images. It is once again seen 
that the performance of the ATI/DPCA algorithms depend strongly on the thresholds, with performance 
that varies across thresholds from image to image. On the other hand, the proposed approach is able to 
detect the targets with high fidelity regardless of the scene/image and does not require tuning of thresholds 
for detection. 
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Fig. 5. This figure compares the output of the proposed algorithm as a function of magnitude and phase for a scene of size 375m by 1200m 
and coherent processing interval of 0.5s. The Bayes SAR algorithm takes the original SAR images in (a) and (b), estimates the nuisance 
parameters such as antenna miscalibrations and clutter covariances, and yields a sparse output for the target component in (c) and (d). In 
contrast, the DPCA and ATI algorithms are very sensitive to the nuisance parameters, which make finding detection thresholds difficult. In 
particular, consider the original interferometric phase image shown in (b). It can be seen that without proper calibration between antennas, 
there is strong spatially-varying antenna gain pattern that makes cancellation of clutter difficult. Calibration is generally not a trivial process, 
but to make fair comparisons to the DPCA/ ATI algorithms, calibration in (f) and (g) is done by using the estimated coefficients from 
the Bayes SAR algorithm. In (e) and (f), the outputs of the DPCA algorithm are applied to the original images (all antennas) and the 
calibrated images (all antennas), respectively. It should be noted that even with calibration, the DPCA outputs contain a huge number of false 
detections in high clutter regions. Nevertheless, proper calibration enables detection of moving targets that are not easily detected without 
calibration, as highlighted by the red boxes. Note that the Bayes SAR algorithm provides an output that is sparse, yet does not require tuning 
of thresholds as required by DPCA and/or ATI. 



TABLE VIII 

Radial velocity estimation (m/s) in 2006 Gotcha collection dataset 



(a) Target 1 



(b) Target 2 



Algorithm 


Bias 


MSB 


No. Missed 


Raw 


0.64 


0.94 


1 


Calibrated 


0.71 


1.02 





Bayes SAR 


0.03 


O.IO 





ATI/DPCA* 


-0.04 


0.20 


27 


ATI/DPCA" 


0.10 


0.20 


2 



Algorithm 


Bias 


MSE 


No. Missed 


Raw 


0.47 


0.77 


6 


Calibrated 


0.48 


0.79 





Bayes SAR 


0.19 


0.22 





ATI/DPCA* 


-0.07 


0.43 


30 


ATI/DPCA" 


0.23 


0.28 


3 



2) Target motion models: Figure 8 shows the output of the proposed approach when prior information 
on the location of the targets might be available. For example, in the shown scene, targets are likely to 
be stopped at an intersection. The performance improvement is given for a mission scene that contains 
target in this high probability region. On the other hand, there are no significant performance decreases in 
the reference scene that does not contain targets in the intersection region. This type of processing could 
be extended to a tracking environment, where targets are projected to likely be in a given location within 
the formed SAR image as discussed in Section IV. 

3) Estimation of radial velocity: The dataset used in this section contained a few GPS-truthed vehicles 
from which we can derive (a) the 'true' location of the target within the formed SAR image, and (b) the 
target's radial velocity which is known to be proportional to the measured interferometric phase of the 
target pixels in an along-track system. Figure 9 shows the estimated radial velocities for two targets over 
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Fig. 6. This figure shows detection performance based on the magnitude of the target response with comparisons between the proposed 
Bayes SAR algorithm and displaced phase center array (DPCA) processing. Note that DCPA declares a detection if the relative magnitude 
to the brightest pixel is greater than some threshold. Results are given for two scenes of size 125m x 125m; within each scene, images were 
formed for two sequential 0.5 second intervals. Scene 1 contains strong clutter in the upper left region, while Scene 2 has relatively little 
clutter. The columns of the figure provide from left-to-right: the magnitude of the original image, the estimated target component from the 
proposed algorithm, the probability of the target occupying a particular pixel, the output of DPCA with a relative threshold of 15 dB, and 
the output of DPCA with a relative threshold of 30 dB. It is seen that DPCA has difficulty in canceling the clutter in Scene 1 with either 
threshold. Moreover, in Scene 2 (c-d) DPCA misses detections of the low-magnitude target in the lower right for the 15 dB threshold. In 
both scenes, there are many false alarms at the 30 dB threshold. On the other hand, the proposed algorithm provides a sparse solution that 
detects all of these targets, while simultaneously providing a estimate of the probability of detection rather than an indicator output. 



18 seconds at 0.25 second increments. We compare the estimation of radial velocity from the output of 
the Bayes SAR algorithm, from the raw images, from the calibrated images, and from two DPCA/ATI 
joint algorithms with phase/magnitude thresholds of (25 deg, 15 dB) and (25 deg, 30 dB) respectively. 
For fair comparisons, the DPCA/ATI thresholds are applied to the calibrated imagery, though this is a 
non-trivial step in general. Numerical results are summarized in Table VIII. It is seen that the Bayes SAR 
algorithm outperforms the others in terms of MSB for both targets. Moreover, the Bayes SAR algorithm 
never misses a target detection in this dataset, which is not the case for the DPCA/ATI algorithms. 

VII. Discussion and future work 

Recent work [20]-[22] has shown that it is possible to successfully decompose natural high-dimensional 
signals/images into low-rank and sparse components in the presence of noise, leading to the so-called 
robust principal component analysis algorithms. [23] introduced a Bayesian formulation of the problem 
that built on the success of these algorithms with the additional benefits of (a) robustness to unknown 
densely distributed noise with noise statistics that can be inferred from the data, (b) convergence speeds 
in real applications of the mean solution that are similar to those of the optimization-based procedures, 
and (c) characterization of the uncertainty (i.e., estimates of the posterior distribution) that could lead 



19 



(Original) 



Z M 



ATI (25 deg) 



ATI (25 deg) 
DPCA (15 dB) 



ATI (25 deg) 
DPCA (30 dB) 



(b) 



(0 



















'v,'t...',. 

















* ■» J 


m 







(d) ' ■ . , 

Fig. 7. This figure shows detection performance based on the phase of the target response with comparisons between the proposed algorithm, 
along-track interferometry (ATI) and a mixture algorithm between ATI/DPCA. Results are given for the same two scenes in Figure 6. In all 
cases, we show results for calibrated imagery where are given by the output of the Bayes SAR algorithm, though this step is not trivial. 
The columns of the figure provide from left-to-right: the phase of the image without thresholding, the estimated target phase component 
from the proposed algorithm, the output of ATI with a threshold of 25 degrees, the output of ATI/DPCA with (25 deg, 15 dB) thresholds, 
and the output of ATI/DPCA with (25 deg, 30 dB) thresholds. In contrast to Figure 6, the contributions from the strong clutter are not very 
strong, though there are still numerous false alarms in the ATI and ATI/DPCA outputs. It is seen that the ATI/DPCA combination with 15 
dB magnitude threshold over-sparsifies the solution, missing targets in (b), (c), and (d). On the other hand, the ATI/DPCA combination with 
30 dB magnitude threshold detects these targets, but also includes false alarms in (a) and (b). 




(b) 
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With TMM 
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(d) 
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Fig. 8. This figure compares the performance of our proposed method with and without priors on target signature locations. In this scene, 
targets are likely to be stopped at an intersection as shown by the region in (a). A mission image containing targets is shown in (b) and a 
reference image without targets is shown in (d). The estimated target probabilities are shown in (c) for the mission scene where inference 
was done both with/without a target motion model (TMM). It can be seen that by including the prior information, we are able to detect 
stationary targets that cannot be detected from standard SAR moving target indication algorithms. The estimated target probabilities in the 
reference scene are shown in (e), showing little performance differences when prior information is included in the inference. 
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(b) Target 2 

Fig. 9. This figure plots tlie estimated radial velocities (m/s) for two targets from measured SAR imagery over 18 seconds at 0.25 second 
increments. Radial velocity, which is proportional to the interferometric phase of the pixels from multiple antennas in an along-track SAR 
system, is estimated by computing the average phase of pixels within a region specified by the GPS-given target state (position, velocity). 
We compare the estimation of radial velocity from the output of the Bayes SAR algorithm, from the raw images, from the calibrated images 
(i.e, using the estimated calibration coefficients), and from two DPCA/ATI joint algorithms with phase/magnitude thresholds of (25 deg, 15 
dB) and (25 deg, 30 dB) respectively. For best comparisons, the DPCA/ATI thresholds are applied to the calibrated imagery, though this is 
a non-trivial step in general. The black line provides the GPS provided radial velocities. Numerical results are summarized in Table VIII. It 
is seen that the Bayes SAR algorithm outperforms the others in terms of MSB for both targets. Moreover, the Bayes SAR algorithm never 
misses a target detection in this dataset, which is not the case for the DPCA/ATI algorithms. 

to improvements in subsequent inference. Moreover, the Bayesian formulation is shown to be capable 
of generalization to cases where additional information is available, e.g. spatial/Markov dependencies. 
Future work will include the development of algorithms that exploit the use of a posterior distribution for 
improved performance in a signal processing task, e.g. detection, tracking or classification. In particular, 
we are interested in using algorithms for simultaneously detecting and estimating targets over a sparse 
scene with resource constraints , as well determining the fundamental performance limits of a SAR target 
tracking system. Furthermore, we would also like to consider other generalizations to the SAR image 
model, such as complex target maneuvers, multiple target classes, and explicit tracking of the target 
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phase. 

Appendix A 
Target signature prediction 

In some applications, such as target tracking or sequential detection, we may have access to an estimate 
of the kinematic state of the target(s) of interest, such as position, velocity and acceleration. This may 
be useful for predicting the location of the target at sequential frames. For simplicity, consider a single 
target whose state (r(r), r'(r)) is known with standard errors {ar, o-,-), where r denotes the slow-time (i.e., 
time of the radar pulse). In standard SAR image formation, moving targets tend to appear displaced and 
defocused in as described in the Uterature by Fienup [1] and Jao [7]. Moreover, Jao shows that given the 
radar trajectory (q, q) and the target trajectory (r, r), one can predict the location of the target signature 
within the image p by solving a system of equations that equate Doppler shifts and ranges, respectively, 
at each pulse: 

^[\\p-Q{r)\\2-\\r{r)-q{T)\\,l^^.^0 (31) 

\\p*-q(r)\\,^\\r(T)-q(r)\\,, (32) 

which can be reduced to the simpler system of equations: 

q{r) ■ [p* - q{T)] = [r(r) - g(r)] ■ [r(r) - q(r)] (33) 
\\p*-q{r)\\, = \\r{r)-q{r)\\, (34) 

In practice, the target state (r, r) is unknown or known with some uncertainty. In the latter case, we can 
predict the probable locations of the target signature by one of several methods, including: 

• Monte Carlo estimation of the target signature locations. 

• Gaussian approximation using linearization or the unscented transformation. 

• Analytical approximation as proposed by Newstadt et al. [26]. 

A. Notation 

Following the derivation of Jao [7], we will assume the following notation: 

• r(r) = {rx-,ry,rz) is the position of a point scatterer. 

• q(r) = (g^., qy.,qz) is the position of the radar platform. 

• r(r) = {rx,ry,rz) is the velocity of a point scatterer. 

• ^{'t) = (Qx, QyAz) is the true position of the platform. 

• P = (Pa;;P2/;Pz) is a pixcl location within the image. 

• T represents the slow-time (i.e., pulse of the radar sample). 

B. Deterministic solution 

In the deterministic case, where r, r, q, and q are all known, we can find the pixel p* where the target 
signature will be focused at time r by solving equations (25) and (26). In particular, if we assume that 
z-coordinate is given by a function 

Pz^h{p^,Py), (35) 

then we can give explicit expressions for {px,Py:Pz) in some cases of h. We will focus on the simple 
case where h{px,Py) = zq (i.e, constant elevation), though this can be easily extended to other cases (for 
example, with a depth elevation map). 
To solve the system of equations, let 

a(r) = ||r(r)-q(r)||^ (36) 
/3(r) = q(r) • r(r) - r(r) • (r(r) - q(r)) (37) 
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Then we have 

a(r) = ||p*-q(r)||? 

= {pI - Qx{T)f + {pI - qy{T)Y + {zo - qz{T)y 
and re-arranging equation (25) we have 

/3(r) = q(r) ■ p* 

= QxPI + qyP*y + QzZo 

For this derivation, assume that 7^ 0^. Therefore, solving for p*, we get: 



(38) 



(39) 



P{T)-q,Zo\ [_%] * (40) 

I • ) Pv 



Px = 



Qx J \ Qx/ 
= 70 + 7iPy 

Plugging into equation (38) we get: 

a(r) = (70 + 7iP; - g.(r))' + {pI - Qyir))' + {zo - g.(r))^ (41) 
which can be re-arranged as 

<^yiP;T + byP; + CyiplT = (42) 

where 

«y= (1 + 7?) 

^2/ = 2 (7071 - 7i?x(r) - ?2/(r)) 

= (70 - Qxir)f + {qy{T)f + {zq - qz{r)f - Q;(r) 
= (7o + ^0) - 2 (laQxir) + Zoq,(r)) + Mr)\\l - \\r(r) - q(T 
Then p* is given by the solution of the quadratic equation: 



(43) 



— by ± ^./bl — ittyCy 

and pI is given by equation (40). This solution suggests that the target energy will generally actually 
appear at two locations. However, in most cases only one of these locations will be in the formed SAR 
image. Thus, we generally choose the solution that is closest to the scene center (0,0). 

Finally, we note that equations (40) and (44) provide the pixel location containing the target energy 
at a single pulse time, r. Generally, images are formed by integrating pulses over a coherent processing 
interval (CPI) containing multiple times re [To,Ti]. 

^By assumption, the radar has non-zero velocity in the xy-plane. Thus, if qx = 0, then this derivation should be vaUd if we switch the x 
and y indices. 
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C. Uncertainty model 

It is unlikely that we will have perfect information regarding the target state r(r) and /"(r) at all times 
r. On the other hand, there are special cases where we might have some information about these states 
that include 

• A tracking environment, where we estimate target position and velocities with standard errors, cr^ 
and (7r, respectively. 

• A persistently monitored scene, where we have knowledge of traffic patterns or road systems. In such 
a case, we have prior knowledge of likely target states. 

In either case, we have a characterization of likely behavior of the target kinematic state. We can represent 
this knowledge in many ways that could include 

• A linear kinematic model, where 

r(r) = ro + VT + 

vr^N (/i,„ all) 
a^N{fi,,all) 

• A random kinematic model, where at each time r 

r{r)r.N{^r.{r),[ar{r)ri) 

Note that both models assume that the position and velocity vectors are Gaussian distributed. However, 
the first model is characterized by only 6 random variables (2 each for position, velocity, and acceleration) 
regardless of the number of pulses. The second model, on the other hand, assumes 4 random variables 
for each pulse r. In fact, the first model can be seen as a specialization of the second model for specific 
structures for the mean and variance parameters. 

In this document, the choice of target kinematic model depends on the inference method which we will 
use to derive the distribution of the target locations. In Monte Carlo sampling, the choice of model is of 
relatively insignificant computational burden as compared to the generation of the Monte Carlo samples. 
On the other hand, in the analytical approximation methods, the choice of target kinematic model is of 
great importance. 

The goal of this section is to provide a prediction model for the locations of targets within a SAR 
image given a target kinematic model. In particular, we define this model through the distribution of pixel 
locations: 

f{Px,Py) (47) 

which is assumed to have support on M^. Generally, images are formed on a discrete grid so that we 
should really consider a discrete distribution. However, for simplicity we consider a continuous domain 
in this section. 

Moreover, we have to be careful how we define the distribution of pixel locations for SAR images 
formed by integrating multiple radar pulses. Consider a probability distribution function (PDF) for the 
target location at pulse Tj given by f{Px{Ti),Py{Ti))- We define the target distribution of interest as: 

1 ^ 

f{Px,Py) = ^^f{Px{ri),Py{Ti)) (48) 

i=l 



This is equivalent to the distribution of the target occupying location {px,Py) at any of T integrated 
pulses. We could consider a richer description by solving for the joint distribution on {Pxi'T'i) , Py{Ti)}J^y 
However, this will be generally very high-dimensional and might not provide any additional benefit over 
the distribution given by equation (48) be useful for the purposes described in the paper. 
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Finally, since we will solve for the distributions f{px{Ti),Py{Ti)) independently for each r^, we will 

consider the random kinematic model only. In the naive situation where the distributions of target locations 
are independent over time, this will provide solutions that can be approximated analytically with just a 
few minor assumptions. 

D. Monte Carlo prediction 

The most straightforward way to approximate the distribution in equation (48) is to use Monte Carlo 
sampling from the linear/random target kinematic models, followed by projection of those target states 
into the image domain using equations (40) and (44). A Monte Carlo representation is subsequently given 
by the average number of samples occupying any pixel. Note that since the Monte Carlo representation 
will contain discrete samples, we will end up with a discrete probability mass function (PMF) rather than 
a PDF. 

E. Gaussian approximation 

Rather than using a potentially high-dimensional PMF or PDF representation of the target location at 
time Tj, we could consider a Gaussian approximation that represents the probability distribution with just 
two parameters: the mean Hi e and the covariance Sj e M^^^ for each slow-time. Then the PDF is 
given by a Gaussian mixture model of form: 



where (pcAf {x; /i, S) is the multivariate normal distribution PDF of x with mean fi and covariance E. 
To find the means /ij and covariances Sj, we could consider linearization of the solution to equations 

(40) and (44) around a particular state ^(tj) = {r{Ti),f{Ti)}. This is akin to the approximation made by 
the Extended Kalman Filter, where the state is Gaussian but the observations are non-linear. Let us define 



Since this is a linear function of a Gaussian distributed vector in ^(rj), we know that the pixels p{Ti) ^ 
G{^{Ti)) are distributed as 




(49) 




(52) 



where S^(ri) is given by 




Sr(ri) 
S,(t,) 



(53) 
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F. Analytical approximation 

It is also possible to get a closer approximation than the linearization example provided above by doing 
some analytical derivations. In particular, let us assume that the radar platform moves in the x-direction 
so that % = qz = ^ and > 0. In the general case, this derivation would hold for a transformed set of 
coordinates {p'^^Py), though we won't go into that derivation here. In the former case, equations (40) and 
(44) reduce to: 

^ f3{T) ^ QxTx - {tx - qx)vx - {ry - qy)Vy ^^^^ 

qx qx 

P*y^qy± V oi{t) - {p^ - q^Y - {zq - q^Y (55) 
From equation (54), we see that 

fipl\r.,ry)^N{ii,a^) (56) 

where 

l^vx\'^x l^vy^y ~l~ qxl^vx ~l~ qyl^vy^ (57) 

qx 

" = Wf 

where we have assumed that S^. = all. Note that cr^ is a function of the position r. However, since 
||r|| in general, we make a zero-th order approximation here so that 

A {li-rx — qx Y + { l^ry — qyY 
{QxY 

In this case, we see that we can find 



^2 ^ ^2 ^ Kt-rx HXJ > ,^ry Hyj ^2 (59^ 



00 00 



f{Px)^ J J f{Px\rx,ry)f{r^,ry)dr:cdry 

(60) 

4>CAf{Px; l^{rx, Ty), (Jo)4>m{rx; l^rx, (^r)<^CAf{ry; Hry, a1)dr^dry 



—00 —00 
00 00 




where we have assumed that = a^I. Since ^{rx-, Vy) is linear in and Ty from equation (57), we can 
analytically solve this integral to see that 



where 



p:~iV(/.,,,aM (61) 



l^px — . {qxfJ'vx ~l~ qyl^vy ~l~ qxf-^rx l^vxl^rx l^vyl^y^ (62) 

qx 

2 

[{qx - I^Vxf + l^ly] + (63) 



a, 



{qxY 

Since p* in equation (55) is non-linear, we make one more assumption with a first-order linearization 
around Note note that given p*, equation (55) only depends on the target state through r (and 

not on the velocity r). Define s{r,Px) to be the value of equation (55) given state r and pixel location 
Px- Then we approximate p* as 

Pl{r)\pl ^ s{tir,P*x) + (ys{fir,P*x)l^^y {r - ^r) (64) 
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Finally, we note that given p^, Py is Gaussian distributed with distribution: 

Py\Px^ N{^ipy,aly) (65) 
flpy = s{fir,Pl) (66) 

Note that both and | have Gaussian distributions that can be described by a mean and covariance 
term. In contrast to Section A-E, the distribution is not jointly Gaussian because the mean and covariance 
of p*y I p* depend on p*. Nevertheless, one can easily evaluate this PDF at any pixel (px^Py)- Over short 
CPIs, both approximations will probably lead to similar results. 



Appendix B 
Inference Details 

In the hierarchical model proposed in Section III, the distribution of hyper-parameters at the base 
layer are generally chosen to be conjugate to the distributions at the next layer. This allows for ef- 
ficient approximation methods for the posterior distribution in the sense that we can sample exactly 
from these distributions. In particular, we use a Markov Chain Monte Carlo (MCMC) algorithm in 
the form of a Gibbs sampler to iteratively estimate the full joint posterior. In MCMC, this distribu- 
tion is approximated by drawing samples iteratively from the conditional distribution of each (random) 
model variable given the most recent estimate of the rest of the variables (which we denote by — ). Let 
= |S, X, G, M, A*^, A^, H, C, T/} represent a current estimate of all of the model variables where 
rj represents the set of all hyper-parameters. Given measurements /, the inference algorithm is given 
in Figure 2. Note that MCMC algorithms require a bum-in period after the Markov chain has become 
stable, where the duration of bum-in period depends on the problem. After this point, we collect N samples 
samples that represent the full joint distribution. 

The sampling details are provided for each of the steps in Figure 2 individually. 



A. Basic Decomposition 

Given the parameters in Tables IV and V, we arrive at one of the primary benefits of using Gibbs Sam- 
pling for inference: namely that we can independently sample across pixels and frames. In distributional 
form, we have 

/(B,X,G,M,A^,A^V,-) (68) 

The conditional independence among pixels and frames given the nuisance parameters allows us to easily 
parallelize the sampling procedure over the largest dimensions of the state. Moreover, we can extend the 
parallelization to samphng independently over passes by separating the sampling of equation (68) into 
two Gibbs steps from the densities: 

f{hf,x%,^^rn^l^^5f^t§\I.-) (69) 

— JVlf I-* ) '^f,l:N-> ''''f,l:N-> "f,l:N ' ) 

i 

f{9%.N.5f^'\l.-) (70) 
^f{5f^\l.-)X[f{g^]\I.5f^\-) 
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It should be noted that each of these distributions have expUcit forms as either multivariate Gaussian or 
Bernoulli distributed. For example consider the distribution 

= /(xg,rng|5j^'(-),rg,-) (71) 

where we define 



■ /hS)-Sf^g^, (72) 
where ./ is the point- wise division operator. In particular, we know that 



r^i-l^^+x^ + S^i^^rr^ + v^, (73) 

Thus, given 5'f-'^\ we have r^} as a linear combination of Gaussian random variables. Thus, by standard 

( ) ( ) 

conditional distributions of a Gaussian random vector, we know that the distribution of bj' and xj- must 
also be Gaussian. Table IX gives the means and covariances of the Gaussian distributions in equations 
(69) and (70). Moreover, the distributions of the indicator functions are easily found by Bayes rule. For 
example, define the quantity 

- f&^^Tf - d,hf},Sf^'\g_fj,bf). (74) 
Then it is easily seen that for the target indicators we have by Bayes rule 

zJjS^/"^ = 1) 

d=m (75) 



^.(p) _ 1 1 r ^ _ ^"f'i 



Note that za is just an evaluation of the Normal PDF so that it can be simply calculated. Table X provides 
explicit values of Zd for both the target and glint indicators, where (l)cj\r{x; /u, S) is the multivariate normal 
PDF with mean /u and covariance S. 



B. Calibration coefficients 

For this thesis, we assumed that pixels within a subset Zg C {1, 2, . . . , P} share the same calibration 
constant so that 

^S, = ZkjAg), Vp e Zg. (76) 

with Zkj,i{g) ~ CM (l, (cr^)^). In our formulation (and dropping the {k,f,i) indices for simplicity) we 
have measurements of the form 

= z{g){&^ + s^P^ + Vj9 e Zg. {11) 

Define = + s^p) which is a known quantity in our Gibbs sampling inference step. Moreover, we 
assume that \y'^P^\ » \v^p^\ so that for any p E Zg we have 

= z{g) 

^ z^P^ + E[z{g)]v^P^ (78) 
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Note that given y^^^ (as in the Gibbs sampling step), we have the situation where z{g), and v^^^ are 
all Gaussian distributed random variables. Thus, the conditional distribution of z{g) is also Gaussian with: 

z{g)\yr.CJ\f{li,{g),a^,(g)) 



1 + 



y) 



V\2 



+ y"y[(y 



H\2 



(79) 
(80) 

(81) 



where y = {y^^^}^^^ and i = {i^^^}^^^^ ■ Note that when (cr^)^ is large, then maximum likelihood 



inference in this case yields the least-squares solution for z{g). 



C. Object class assignment 

In this model, we assume that each pixel can be assigned to one of J possible classes. We assume that 
the number J is known a priori and do not consider the details involved in the merging or splitting of 
object classes here. More detailed models (such as the so-called Indian Buffet processes) can also estimate 
the number of classes directly from the data. 

In this model, inference on class assignment is straightforward given the distributions (i.e., covariance 
matrices) for each class. Define the matrices 



^2 



'ip) 



e C 



FxK 



X 



X 



ip) 
1,1 

(p) 

1,2 



e C 



FNxK 



(82) 



Then the probability that pixel p belongs to class j is given by: 

^(p) ^ Pr(pixel p has class j) = exp{TB + Tx + Qj} 
where qj is the prior probability of class j and 

Ts = -trace {[T^ {j)]-'b^^\b^^Y) - ^ log |r^(i)| - KF\og{2n) 

Tx = -trace {[V'' {j)]-'x^p\x^pY) - — log |r^(j)| - i^FiVlog(27r) 
Then the class assignment to pixel p is the single location in c^^ with value equal to one, where 



Multinomial(l;w(^)) 



(83) 

(84) 
(85) 

(86) 



Note that we can improve upon this model by allowing the probabiUties for pixel p to vary spatially 
(i.e., pixels are likely to share the same class with neighboring pixels). One simple way to include this 
information is to let 

^(P) ^ ^(P) ^ ^(P) (87) 

where g^^^^j^ is some filter for averaging nearby pixels and * is the convolution operator (assumed to be 
supported on the same set of pixels as 'u^'^^'). Then we draw 

c(p) - Multinomial(l;m(^')) (88) 



D. Hyper-parameters 

In this model, we have three types of hyper-parameters that need to be estimated: covariance matrices 
(or variances), indicator probabilities, and object class probabilities. In all cases, the distribution of these 
parameters depend on test statistics of much smaller dimension that P. 
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1 ) Covariance matrix inference: We model the covariance matrices for the Normal distributions in two 

ways: (1) for the stationary components (background, speckle, and glints), we model the covariance matrix 
as a random variable; and (2) for the other components (targets, additive noise, calibration coefficients), 
we assume independence among the antennas. In particular, consider a random vector of K elements, w_, 
with 

wr^CN (0, a'^Tp) (89) 
(7^ ~ InvGamma(aCT, 6o-) (90) 

Then in the stationary case, we have 

Tp ~ InvWishart (ar((l - p)Ikxk + P^rIk), i^r) (91) 
p ~ Beta(ap, bp) (92) 

and in the independent case, we have 

Tp = Ikxk (93) 

First consider the case where T p is a random variable. Assume that we have n independent samples of 
ui, which we refer to as 1^ = vec {w}. Then, we consider a Gibbs sampling procedure: 

Sample ~ f{(y'^\W, Tp, p) (94) 
Sample ~/(r„p|W,(72) (95) 

Let T = 1/(7^ ~ Gamma(ao-, ha). Then 

f{T\W,Tp,p)^f{W\T,Tp)f{T) 



oc 



exp |-^trace(r;i"Vyvr^)} 



r(a,) 



r"" ^exp{-6o-T} 



(96) 



cx: 



-r'^-^xpl-ftV} 



where 

a! — + 



2 

tYBjceiV-^WW^) ^^^^ 

b' = b, + ^-^ ^ 

This demonstrates that in this situation, o"^ has an Inverse-Gamma distribution with parameters a' and b'. 
Note that in the case where Tp = Ikxk, then the posterior parameters are given by 

a' = a, + - (98) 

Thus, in the Gibbs sampling procedure, the variance parameter a'^ is Inverse-Gamma distributed whether 
or not T p is modeled as a random variable. Table XI provides the posterior Inverse Gamma distribution 
parameters for the variance parameters in our model, where vec {■} refers to the vectorization operator. 

For the background, speckle, and glint components, we also need to sample the coherence parameter p 
and correlation matrix T p. Let Vl^ = W/cr be our observed measurements given as given by equation 
(95). Define fx — n/a. Then we have 

T^l {Tp,a\p)r.CN{ix,Tp) 

Tp ~ InvWishart ([pl/f 1^ + (1 - p)lKy.K\{y -K-l),v) (100) 
p ~ Beta(ap, bp) 
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Note that this is in the form of the Multivariate-Normal-Inverse-Wishart conjugate distribution given p. 
This leads to the well known posterior parameters: 

H 

r ■ ■ 



{W, a\ p) - InvWishart(A^ar + (^^ - A) [v^ - A) , i^r + n) (101) 



m=l 



where Ap = + (1 ~ p)Ikxk- Ideally, we would like to sample both Tp and p jointly. Even though 

we can simply sample from equation (101), the same is not true for the density 

p\{w,a^p^ (102) 

which is required in order to jointly sample these parameters. Fortunately, we know that 

/(r„ p\W, a') oc f(Tp\W, a^ p)f(p) (103) 

which is easily evaluated since we have closed form functions for both of these densities. Thus, we can 
use Metropolis-Hastings to sample p and Tp. 

2) Indicator probabilities: In the basic model where the indicator Beta distribution parameters do not 
depend spatially or temporally, then the posterior indicator probabilities for 

6 ~ Bemoulli(7r) (104) 

TT ~ Beta(a7r, b^^) (105) 

are given by 

7r|5~Beta(a^ + 5,6^ + (l-5)) (106) 

Note that we can modify and as in Section IV-A. However, the posterior inference for the probabilities 
is identical by replacing and b.„ by their spatiotemporally varying version. 

3) Object class probabilities: We use a Multinomial-Dirichlet conjugate pair to determine object class 
assignments, where the class probabilities q have a prior Dirichlet distribution with cj = 1/ J for j ~ 
1, 2, . . . , J. Then, after observing the class assignments, we can calculate the number of pixels in any 
class 

Nj^\Qj\^\{p:c^^ ^j}\ (107) 



Then the posterior distribution for the class probabiUties is given by 

'i=i 

where [N]j = Nj. 



g\ {^ilLi ~ Gamma {q + N, l) / J (108) 
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